function   res = momentsClosedForm(model,setupFilter)

% For efficient starting values when computing unconditional second moments
%global Var_z2old 
%if ~exist('Var_z2old','var')
    Var_z2old = [];
%end
% Construct model structure to match the model output
modelGMM = doModelforGMM(model);

% The moments in the normal distribution
ne         = size(model.eta,2);       %number of shocks
vectorMom3 = zeros(1,ne);
vectorMom4 = 3*ones(1,ne);

% Computing moments
maxAuto_lag = setupFilter.CSregress_m;
mom2nd      = UnconditionalMomentsProj_2nd_Lyap(modelGMM.g0,modelGMM.gx,modelGMM.gxx,...
    modelGMM.h0,modelGMM.hx,modelGMM.hxx,modelGMM.eta,vectorMom3,vectorMom4,maxAuto_lag,Var_z2old);
%Var_z2old = mom2nd.Var_z;

if model.appOrder <= 2
    % hours and wages are in deviation from steady when matched to the data
    pos_l      = find(strcmp(modelGMM.labely,'$l_t$'));
    pos_w      = find(strcmp(modelGMM.labely,'$w_t$'));
    mom2nd.E_y(pos_l,1) = mom2nd.E_y(pos_l,1)-modelGMM.gSS(pos_l,1);
    mom2nd.E_y(pos_w,1) = mom2nd.E_y(pos_w,1)-modelGMM.gSS(pos_w,1);
   
    % Transforming moments from log(c_t) to consumption growth, log(c_t)-log(c_t-1) + log(myz_t)
    pos_c           = find(strcmp(modelGMM.labely,'$c_t$'));
    pos_c_ba1       = find(strcmp(modelGMM.labelx,'$c_{t-1}$'));
    pos_myz         = find(strcmp(modelGMM.labelx,'$\mu _{z,t}$'));
    mom2nd.E_y(pos_c,1)  = 4*log(modelGMM.params.MUZss);
    % var(log(c_t)-log(c_t-1) + log(myz_t)) = var(log(c_t))+var(log(c_t-1)) + log(myz_t) 
    %                                        -2Cov(log(c_t),log(c_t-1)) + 2Cov(log(c_t),log(myz_t)
    %                                        -2Cov(log(c_t-1),log(myz_t))
    Cov_y_z  = mom2nd.C*mom2nd.Var_z;
    varC     = mom2nd.Var_y(pos_c,pos_c);
    varC_ba1 = mom2nd.Var_x(pos_c_ba1,pos_c_ba1);
    varMUZ   = mom2nd.Var_x(pos_myz,pos_myz);
    % recall: cov(y,x) = cov(y,xf + xs) = cov(y,xf) cov(y,xs)
    CovC_C_ba1 = Cov_y_z(pos_c,pos_c_ba1)+Cov_y_z(pos_c,modelGMM.nx+pos_c_ba1);
    CovC_MUZ   = Cov_y_z(pos_c,pos_myz);
    CovC_ba1_MUZ = mom2nd.Var_x(pos_c_ba1,pos_myz);
    varDc = varC + varC_ba1 + varMUZ - 2*CovC_C_ba1 + 2*CovC_MUZ - 2*CovC_ba1_MUZ;
    mom2nd.Var_y(pos_c,pos_c)         = 4^2*varDc;
   
    % recall Var(y) = E[y^2]-E[y]E[y'], so E[y^2] = Var(y)+E[y]E[y']
    tmp               = diag(mom2nd.Var_y+mom2nd.E_y*mom2nd.E_y');
    EyEy2             = [mom2nd.E_y(1:setupFilter.nyMom,1);tmp(1:setupFilter.nyMom)];
    
    % Loadings for yields in the extended state space system (with pruning)
    D_yields         = 4*[model.byx  model.byx  model.Byxx];
    
    % Loadings for the slope in the CS regressions in the extended state space system (with pruning)
    D_slopeCS        =  D_yields(setupFilter.CSregress_m+1:end,:)-D_yields(setupFilter.CSregress_m,:);
    
    % Getting Cov(yields^k,slope^k)
    Cov_y_k_slope_k = D_yields(setupFilter.CSregress_m+1:end,:)*mom2nd.Var_z*D_slopeCS';

    % Getting Cov(yields^(k-m),slope^k). Recall Cov_z(:,:,m) = Cov(z_t+m,z_t)
    Cov_y_kminusM_slope_k  = D_yields(1:end-setupFilter.CSregress_m,:)*...
                             mom2nd.Cov_z(:,:,setupFilter.CSregress_m)*D_slopeCS';

    % Getting Var(slope^k) 
    Var_slope              = D_slopeCS*mom2nd.Var_z*D_slopeCS';                         
    
    % First and second moments of all yields
    meanYields             = 4*model.by0 + D_yields*mom2nd.E_z;
    stdYields              = sqrt(diag(D_yields*mom2nd.Var_z*D_yields'));
    
    % First and second moments of shrinkage moments used in estimation
    meanForEstim           = mom2nd.E_y;
    stdForEstim            = sqrt(diag(mom2nd.Var_y));
elseif setupFilter.orderApp == 3
    error('momentClosedForm.m: no closed-form moments for third order approx')
end

Cov_excessReturn_slope = Cov_y_kminusM_slope_k-Cov_y_k_slope_k;
idx                    = setupFilter.CSregress_m+1:size(D_yields,1);
scaling                = setupFilter.CSregress_m./(idx-setupFilter.CSregress_m);
CS_cov                 = [NaN(setupFilter.CSregress_m,1);scaling'.*diag(Cov_excessReturn_slope)];
CS_loadings            = [NaN(setupFilter.CSregress_m,1);diag(Cov_excessReturn_slope)./(scaling'.*diag(Var_slope))];
%% Output

if setupFilter.SMMwithCSreg == 2
    % Including CS moments and adjusted CS moments
    %risk-adj CS moments most be the empirical moments and therefore added later
    CS_var           = [NaN(setupFilter.CSregress_m,1);scaling'.^2.*diag(Var_slope)];
    res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS);...
                        CS_var(setupFilter.maturites_CS)]';
elseif setupFilter.SMMwithCSreg == 1
    %res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS)]';
    CS_var           = [NaN(setupFilter.CSregress_m,1);scaling'.^2.*diag(Var_slope)];
    res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS);CS_var(setupFilter.maturites_CS)]';
else
    res.modelMoments = EyEy2';
end
% Removing first moment for mean of hours and wages
res.modelMoments    = res.modelMoments(3:end);

res.CS_cov           = CS_cov;
res.CS_loadings      = CS_loadings(setupFilter.CSregress_m:setupFilter.CSregress_m:end)';
res.meanYields       = meanYields;
res.stdYields        = stdYields;
res.meanForEstim     = meanForEstim;
res.stdForEstim      = stdForEstim;
scalingAll           = [NaN(1,setupFilter.CSregress_m) scaling];
res.scaling          = scalingAll(setupFilter.maturites_CS);
    



end